Cooling 



of mechanical motion with a two level system: the high temperature regime 
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We analyze cooling of a nano-mechanical resonator coupled to a dissipative solid state two level 
system focusing on the regime of high initial temperatures. We derive an effective Fokker-Planck 
equation for the mechanical mode which accounts for saturation and other non-linear effects and 
allows us to study the cooling dynamics of the resonator mode for arbitrary occupation numbers. 
We find a degrading of the cooling rates and eventually a breakdown of cooling at very high initial 
temperatures and discuss the dependence of these effects on various system parameters. Our results 
apply to most solid state systems which have been proposed for cooling a mechanical resonator 
including quantum dots, superconducting qubits and electronic spin qubits. 
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Over the past years quantum ground state cooling 
schemes for micro- and nano-mechanical resonators have 
attracted considerable attention. Work in this direc- 
tion is motivated by the goal to study quantum effects 
with macroscopic objects [U-0 as well as potential ap- 
plications for quantum information processing and 
nano-scale sensing techniques ■ Apart from passively 
cooling high frequency mechanical modes in a cryogenic 
environment |lOj several active cooling schemes are cur- 
rently explored. So far the most successful techniques 
are based on opto-mechanical cavity cooling ideas jlll - 
Il7j where the mechanical resonator is damped by inter- 
actions with a driven optical fl8l - [24| or microwave (25l[26j 
cavity. Alternatively, mechanical motion can be cooled 
by coupling the resonator to a dissipative two level sys- 
tem (TLS) and various implementations including quan- 
tum dots [27|, [28|, superconducting qubit s I29l -l34j and 
nitrogen-vacany (NV) defects in diamond [35| have been 
suggested in this context. Although in practice more 
challenging such an approach offers the great advantage 
that once the resonator is cooled to the ground state, the 
TLS can be used as a non-linear element to prepare and 
detect non-classical states of the resonator mode. 

For cavity cooling schemes the linearized opto- 
mechanical interactions allow for an exact analytic treat- 
ment of the cooling problem [TH, [ljj which in general 
is not the case for a mechanical resonator coupled to a 
TLS. Therefore, most theoretical studies of TLS-based 
cooling schemes (27rl29l I33rl35| are focused on the regime 
of weak resonator-TLS interactions which is analogous 
to the Lamb-Dicke (LD) regime discussed in the context 
of laser cooling of trapped atoms and ions 36-38]. The 
cooling dynamics is then described by a rate equation for 
the mean resonator occupation number n r (t), 



.(t) = -T(n r (t) - n f ), 



(1) 



where the damping rate T and the final occupation num- 
ber rif can be calculated from the linear response of 
the TLS. In this regime opto-mechanical and TLS-based 
cooling schemes lead to qualitatively similar predictions 
for F and n / which in essence also agree with the known 
results from laser cooling of trapped atoms [36-38]. 



Although the weak coupling approximation is usually 
well justified at low resonator occupation numbers a sim- 
ple rate equation ([T]) cannot be valid for arbitrary tem- 
peratures since the maximum cooling power is naturally 
bound by the decay rate of the TLS. Recently we have 
proposed to use optical pumping of the spin states of a 
NV center in diamond to cool the motion of nano-scale 
cantilevers 35]. This setup can in principle be operated 
at room temperature where for oscillation frequencies 
of u r /2Tr 1 MHz the thermal equilibrium occupation 
number can be as high as Nth ~ 10 6 - Especially in this 
extreme case one can expect that non-linear effects such 
as the saturation of the TLS lead to considerable modi- 
fication of the cooling dynamics. However, also for other 
physical implementations and lower initial temperatures 
of T = 0.1 — IK the resonator energy can still vary over 
several orders of magnitude during the cooling processes 
and the validity of Eq. ([1]) under different experimental 
conditions must be addressed. 

In this work we analyze the cooling dynamics of a me- 
chanical resonator mode coupled to a dissipative TLS, 
focusing in particular on the regime of high initial tem- 
peratures. Based on a semi-classical approximation we 
derive an effective Fokker-Planck equation for the res- 
onator mode which is valid for arbitrary occupation num- 
bers and allows us to study the crossover from the LD 
to the high temperature regime. Already at moderately 
high occupation numbers we find a degrading of the cool- 
ing rates, but in this case the final state is still accurately 
described by the LD theory. However, beyond a certain 
critical value of Nth cooling is highly suppressed and can- 
not compete with environment-induced re-thermalization 
processes. In this regime cooling of a mechanical mode 
with a TLS is no longer possible. We study the crossover 
between the different cooling regimes and discuss the de- 
pendence of the critical occupation number on the rele- 
vant system parameters. 

The paper is structured as follows. In Sec. U we intro- 
duce the basic model for a mechanical resonator coupled 
to a solid state TLS and present in Sec.|H]a brief overview 
of the cooling results in the LD regime. In Sec. IIIII we 
then study the dependence of the cooling rate on the res- 
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onator occupation number and derive in Sec. lIVI an effec- 
tive Fokker-Planck equation to describe the fuli cooling 
dynamics and steady state properties of the system. In 
Sec. [V] we estimate the effects of higher order mechani- 
cal modes and finally in Sec. I VII we summarize the main 
results and conclusions of this work. 



I. MODEL 

We consider a nano-mechanical resonator coupled 
to a dissipative TLS, for example a superconducting 
qubit |2SH34|. a quantum dot [27| or an electronic spin 
qubit j35j. We model the resonator as a single harmonic 
mode with frequency u r and annihilation (creation) oper- 
ator a (a^). The validity of this single mode approxima- 
tion and the influence of higher order vibrational modes 
will be addressed in Sec.|V] The lower and upper state |0) 
and |1) of the TLS are split by a bare transition frequency 
woi and coupled by a classical driving field of frequency 
uid = ^01 + S. The motion of the resonator leads to an 
additional modulation of the TLS splitting and in the 
frame rotating with Ud the total system Hamiltonian is 

(n = i) 

6 fl t A t 

H = --j°z + -^o-x +u) r a l a+ -^{a + a [ )(T z . (2) 

Here at are the usual Pauli operators, a (a^) the anni- 
hilation (creation) operators for a mechanical mode of 
frequency uj r and il is the Rabi frequency of the classi- 
cal driving field. The last term in Eq. ([2]) describes the 
TLS-resonator coupling where A is the frequency shift 
per vibrational quanta. Depending on the specific physi- 
cal implementation the origin of this frequency shift can 
be due to electrostatic, magnetic or deformation poten- 
tial interactions, and for a more detailed derivation of 
Hamiltonian ([2]) we refer the reader to the original pro- 
posal cited above. 

A. Master equation 

The TLS and the resonator mode interact with the 
environment and we model the resulting dissipative dy- 
namics by a master equation 

p=-i[H,p]+£ r (p)+£ 1 (p), (3) 

for the system density operator p. The Liouville operator 
Ct accounts for relaxation and dephasing processes of the 
TLS and is in general given by 

£r(p) = — (NTLS + 1) (2(7-/9(7+ - <T + (T-p - p<T+<T-) 

+— TVtls (2(7+ per _ - (T-cr+p - per- (J + ) 

+ Y {VzPOz - P) ■ 

(4) 



Here the first two lines describe relaxation of the 
TLS towards the equilibrium value (<J z ){t — > 00) = 
-l/(2iV T LS + 1) with a total rate Tf 1 = r ± (2iV TL s + 1)- 
Note that the effective occupation number may differ 
from the thermal equilibrium value iVrLS = (e huol ^ kBT — 
l)- 1 in the presence of experimental imperfections. The 
third line accounts for additional dephasing processes and 
the total dephasing rate of the TLS is T 2 _1 = Tf 1 / 2 + r ||- 
In case the TLS is represented by the spin states of a NV 
center the effective decay rate T± is adjustable by the 
optical spin pumping rate and F||,./Vtls > are related 
to non- ideal optical pumping processes [35j . 

The third term in Eq. ([3]) describes mechanical dissi- 
pation due interactions with a thermal phonon reservoir 
of the support or other intrinsic damping mechanisms. 
We model mechanical dissipation by 

£"y{p) =—{N t h + 1) (2apa) — a)ap — pa? a) 

it (5) 

-iVth (2a) pa - aa) p - paa)) , 

where 7 = u r /Q is the energy damping rate for a me- 
chanical quality factor Q and N th = (^"v/feeT _ ^-1 j s 
the equilibrium phonon occupation number. Throughout 
this work we focus on the high Q and high temperature 
regime where "fN t h — ksT/hQ is independent of the res- 
onator frequency. 

B. Displaced oscillator basis & the Lamb-Dicke 
parameter 

Our goal in the following is to study the effective cool- 
ing dynamics of the resonator mode which results from 
interactions with the dissipative TSL. Before we proceed 
let us briefly remark on the relation between the TLS- 
resonator coupling defined by Eq. ^ and the basic model 
for atom-light interactions, 

Ha l = 0£ ^ e inLn^+a) a+ + e -ir, LD ^+a) a _^ ^ ^ 

discussed in the context of laser cooling of trapped atoms 
and ions [3&14381 ] . Eq. (0 describes the interaction of a 
trapped two level atom with a laser of Rabi frequency 
flL- The displacement operators account for the pho- 
ton recoil whenever the atomic state is flipped. Here 
T]ld = k^ao is the Lamb-Dicke parameter for a laser 
wavevector and ao the size of the groundstate wave- 
function of the trapped atom. It can be re-expressed as 
T]ld = \J ^rccoii/^* where w reC oil ~ Tik 2 L /2m is the recoil 
and u>t the trapping frequency. For a tight confinement 
and low occupation numbers, t\ld\J (a) a) <C 1, Eq. ([6]) 
can be expanded to lowest order in t]ld which, for exam- 
ple, leads to a simplified description of sideband cooling 
schemes for trapped ions [37l 38] . For weak confinement 
or high temperatures this expansion is no longer possi- 
ble. In this case e = w rcco ii/r+ = rf^ D uJt/^± has been 
identified as the relevant expansion parameter which for 
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£ < 1 still allows for a semiclassical treatment of laser 
cooling problems (36j . 

As pointed out, e.g., in Ref. the connection be- 
tween Eq. ([2]) and Eq. (|6|) is most apparent in the dis- 
placed oscillator basis which is defined by the unitary 
(polaron) transformation H = UHW where U = e lS 
and S = —i\/(2uj r )(a' ! — a)a z . In this basis we obtain 



H 



ft / -*-( a t-a) 



(7) 

and see that in analogy to the momentum kick in Eq. ([6]) 
a flip of the TLS state is now accompanied by a displace- 
ment of the resonator state in position space. By com- 
paring Eq. ([7]) and Eq. we find that in the present 
system r) = \/uj r is the equivalent of the LD parameter 
while the quantity X 2 /uj r can be identified with the re- 
coil frequency. Therefore, we can already at this stage 
argue that for the experimentally most relevant regime 
A <C uj ri T±_ an approximate treatment of the system 
dynamics in close analogy to semiclassical laser cooling 
problems should also be applicable for the present sys- 
tem. 



II. THE LAMB-DICKE REGIME 

As a starting point we first consider in this section 
the weak coupling and low temperature regime where 
^y/Nth + 1/2 -C Tj_,uj r is satisfied. As discussed above, 
this limit is equivalent to the LD regime in laser cool- 
ing and has also been analyzed in several previous works 
on cooling schemes for mechanical resonators [2714291 133| - 
l35j |. Here we only present a brief summary of the main 
results in this limit which we will be the basis for further 
discussions in Sec. IHII and Sec. IIV1 



A. Rate equation in the LD regime 

In the LD regime the effect of the resonator-TLS in- 
teraction can be treated as a small perturbation and it 
is convenient to decompose the master equation Q into 
three terms 



P = £tls(/?) + £ r (p) + £\(p), 



(8) 



such that £tls and C r describe the unperturbed evolu- 
tion of the TLS and the resonator respectively and 



(9) 



accounts for the coupling. In Eq. © we have set A<r z — 
&z — (cz)o and omitted a constant force proportional to 
the steady state value (cr z )o- This force does not con- 
tribute to cooling and can be reabsorbed in a shift of the 
resonator equilibrium position. Under the action of £tls 
the TLS relaxes to its steady state p^LS 011 a time scale T\ 



which is fast compared to the timescale associated with 
C\. Therefore, the total system density operator can be 
to a good approximation written as p{t) ~ Ptls ® Prif) 
and the effects of C\ on the dynamics of the resonator 
state pr (t) can be evaluated in second order perturbation 
theory [37[ . For the present system this calculation is de- 
tailed in Ref. [Ill] and as a result we obtain an effective 
master equation which in the frame rotating with uj r can 
be written as 

r 

p r = Cj(p r )+-^(No + 1) (2ap r a) — a)ap r — p r a*a) 



H — — Nq (2a^ p r a — aa) p r — p r aa^ 



(10) 



Here we have introduced a cooling rate F c = S(uj r ) — 
S(—u> r ) and the minimal occupation number Nq — 
S(—u) r )/T c which are determined by the equilibrium fluc- 
tuation spectrum 



A 2 f°° 
S(u) = —ReJ^ dT((a z (T)a z (0)) Q 



(11) 



From the effective master equation (|10p we can now de- 
rive the rate equation (TTJ) for the resonator occupation 
number n r (t) = Tr{a^ap r (t)} with a total damping rate 
r = T c + 7 and rif — tild- Here 



jN th 



+ N , 



(12) 



is the steady state occupation number in the LD limit. 

In the LD regime the cooling dynamics of the resonator 
mode is fully determined by the fluctuation spectrum 
S(±oj r ) which characterizes the ability of the unperturbed 
TLS to absorb or emit energy at the mechanical fre- 
quency ui r . Fig. [T] shows the typical behavior of S(cu) 
together with the energy level diagram which qualita- 
tively explains the relevant cooling and heating processes. 
The resonance at uj = uj ge = \J S 1 + ft 2 corresponds to a 
transition between \g) — > |e) which are the two dressed 
eigenstates of the driven TLS, 

\g) = cos(0/2)|O)-sin(0/2)|l), 
|e) = sin(0/2)|O) +cos(0/2)|l>, 

where 8 — arctan(ft/|<5|). This process leads to cooling 
and in the sideband resolved regime u} ge T2 3> 1 it is most 
effective when oj r ~ uj eg . For a finite steady state popula- 
tion in the state |e) the reverse process ~ S(~u} eg ) leads 
to heating. This occurs for iVrLS^ii > or under very 
strong driving conditions ft 3> S. In the Doppler regime 
TiU) ge < 1 the cooling and heating resonances start to 
overlap and also pure diffusion processes associated from 
the resonance at uj — become important. In this regime 
cooling is less efficient, but as we see below also more ro- 
bust. 

The spectrum S(uj) can be evaluated using the quan- 
tum regression theorem [§9| which allows us to calculate 
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FIG. 1: a) Fluctuation spectrum S(oj) where uj ge — y/8 2 + Q 2 
is the frequency splitting between the two dressed states \g) 
and \e). The solid line shows the result for Q = 0.75o; r , 
Fx = O.lcJr and the ideal case = JVtls = 0. The other 
two curves indicate the effects of a finite thermal population 
TVtls = 0.15 (dashed line) and excess dephasing Ty = 2Fj„ 
(dotted line), b) Level scheme of the coupled resonator- TLS 
under resonance conditions uj r = U) ge , where \n) denote the 
vibrational eigenstates. 



the cooling rate r c and the minimal occupation num- 
ber No for arbitrary system parameters. The general 
expressions are given in App. [S] and in the following we 
summarize the main results. 



B. Final occupation numbers 

In the limit of an isolated resonator mode 7 — > the 
final occupation number is limited by ii^c ~ Nq. For a 
pure decay process, i.e. T\\, Ntls — > 0, we find that N 
is independent of the Rabi-frcqucncy Q and by choosing 
an optimal detuning S = — y/urf + / 4 we obtain 

In the two limiting cases Tj_ -C U) r and T ±_ 3> U) r 
Eq. (fT5|) reproduces the well known results for sideband 
and Doppler cooling respectively. In particular, this 
means that T± < uj t is a minimal requirement to achieve 
ground state cooling. In a solid state environment excess 
dephasing processes T|| > are often non-negligible and 
cause an additional broadening of the TSL transition. 
While in this case the detailed behavior of No is more 
complicated, we find that in essence ground state cool- 
ing requirements are now determined by the condition 
T2U> r > 1 as expected from the qualitative discus- 
sion given above. Finally, finite temperature effects or 



other imperfections can cause a non-vanishing occupa- 
tion number JVtls which introduces an additional limit 
-/Vo > A^tls on the minimal occupation number. More 
over, we find that for weak driving and sideband resolved 
conditions 

A^o « iV TLS + + /V TLS ll + 0(N 2 LS ), (14) 

which leads to a divergent occupation number for 57 — ¥ 
0. This can be understood from the fact that thermally 
activated jumps between the states |0) and |1) cause a 
fluctuating force, which for £1 — > is not compensated 
by a corresponding damping mechanism. This effect does 
usually not appear in atomic laser cooling, but has been 
discussed, ejj., for microwave cooling schemes for polar 
molecules [4fJ. 

In summary we see that ground state cooling can in 
principle be achieved under sideband resolved conditions 
IT j_ , Ty < LJ r , low occupation numbers ./Vtls < 1 and a 
minimal driving strength f2 > -/VtlsTj^. For most imple- 
mentations discussed in the literature these conditions 
can be satisfied. Therefore, we expect that in many ini- 
tial experiments the final occupation number will not 
be limited by No, but rather by the competition be- 
tween cooling and thermalization processes. In this case 
tild — jNth/^c and therefore in the rest of the paper we 
will mainly focus on the cooling rate T c . 

C. Cooling rates 

Let us now study the cooling rate T c which in gen- 
eral is maximized for different parameters than No is 
minimized. For the following discussion we can as- 
sume iVxLS ~ p~fl ] and we introduce a parameter 
e = 2T1/T2 > 1 which characterizes the excess dephasing 
compared to the bare decay rate. We have already ar- 
gued above that in the sideband resolved regime T 2 w r > 1 
cooling is optimized when u> r matches the eigenfrequency 
ujge = y/S 2 + fl 2 of the dressed TLS. Therefore, we can 
set 5 = —Lu r cos(6>) and fl — uj t sin(0) for 6 G [0, n/2]. We 
then obtain 

_ A 2 2 sin(0) sin(20) 

c ~ rT (2 cos 2 (6») + e sin 2 (0) ) (2 sin 2 (0) + e + e cos 2 (0) ) ' 

r4\ 5) 

and similar expressions have been derived in Ref. [33| for 
the case e = 1 and in Ref. (3f| for a three level system. 
In Fig. [2] a) we plot T c for different parameters. The 
maximal cooling rate of T c w 0.44 x A 2 /r^ is achieved 
for e — 1 and w 7r/3, which corresponds to a driving 
strength of fi sa 0.850;,- [33j. 

In the Doppler regime Tiio r -C 1 the transition fre- 
quency of the TLS can no longer be resolved and in 
this case the cooling rate is optimized for a detuning 
6 w -1/(4T 2 ). In Fig. Hb) we plot T c in the Doppler 
regime as a function of fi. We find that compared to the 
sideband limit the maximal cooling rate is reduced by 
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FIG. 2: Cooling rates in the LD regime, a) Cooling rate r c in 
the sideband limit u r T2 1 for different e = 2T1/T2 and cj r — 
\fip~- VTW. b) Cooling rate F c in the Doppler limit cj r T2 <S 
1. For each value of Q. and e the detuning is numerically 
optimized as indicated by the dashed line for e = 1. c) Cooling 
rate T c as a function of T± and SI for the ideal case e — 1. For 
each parameter set the detuning 8 is optimized to maximize 

r c . 



another factor of w r T2 and the maximal cooling rate is 
achieved for a driving strength Q « T^ 1 . Finally, Fig. [5] 
c) shows the cooling rate as a function of T± and n for the 
ideal case e = 1. From this plot we conclude that cool- 
ing rates close to the maximum value of T c ~ 0(A 2 /rjJ 
can be achieved for decay rates up to Tj_ ~ but only 
under strong driving conditions Q ~ max{u r , T^}. 



III. COOLING RATES IN THE HIGH 
TEMPERATURE REGIME 

Our goal is now to extend the previous cooling results 
to the regime of high temperatures where the condition 
X\fN t h "C r_i_ is violated. As a first step we will in 
this section use a simple semi-classical approximation to 
study the dependence of the cooling rate on the mean 
resonator occupation number. A more rigorous discus- 
sion of the full cooling dynamics and the steady state of 
the system is then presented in Sec. IIVI 



A. Semi-classical cooling rates 

Let us change into a frame rotating with uj r and for 
the moment assume that the resonator is prepared in a 
coherent state p r = \at){a\ with \a\ 2 ~ Nth ^> 1. In the 
limit 7 — > the change in the occupation number is then 



given by 

d t (a^a) = -i± (aV*^ - ae~^) (a z (t)). (16) 

We introduce the Bloch vector S = (<j-,cr +7 <j z ) T and 
evaluate its dynamics using a semi-classical approxima- 
tion, e.g. (acrfc) ~ (a)(CTfc). As a result we obtain modified 
Bloch equation which in a matrix notation can be written 

as 

(S) = A0) - iX (e- iUrt a + e iu ^a*) A Z (S) - r ± V z , 

(17) 

where V z = (0,0, 1) T and (A,) u = -(A,) 2a = 1 and 
(A z )ij = otherwise. The first term in Eq. (fl7|) is the 
free TLS evolution where 

/i5-jr iQ/2 \ 
A = 2 -iS- -in/2 , (18) 

\ m -in 2 -i / 

and has eigenvalues of 0(T{~ ). The second term in 
Eq. (fTT)) describes an additional driving field ~ Xa 
caused by the resonator motion. To proceed we use the 
continued fraction method developed by Stenholm and 
Lamb [42], |43| and write the Bloch vector in a Floquet 
representation 

00 

(S)(t) = S n {t)e~ inUrt . (19) 

n— — 00 

Here S n (t) are slowly varying coefficients which relax to 
their steady state values Sft on a timescale Tj . For XT\ <C 
1 the resonator amplitude does not considerably change 
on this timescale and in Eq. we can treat a as a 
constant parameter. Then, by inserting Eq. (|19l) into 
Eq. (fTTj) we can express the steady state values of S z \ (a) 
in terms of a matrix continued fraction 

-iaXA z S r l - 1 + (A+iuj r nl)S%-ia*XA z SZ +1 = S nQ T ± V^ 

(20) 

which can be efficiently evaluated numerically. After in- 
serting the steady state values back into Eq. (fTB) and 
taking the average over a few oscillation periods we ob- 
tain 

d t {<Ja) ~ -i± (a*S+o(«) " ^J(a)) = -r c (a)|a| 2 , 

where we have defined the amplitude dependent cooling 
rate T c (a) := iX(Sl /a — S~q/oi*) which depends on \a\ 
only. Note that in deriving Eq. (|2"Tj) we have assumed 
that A -C Tj_, but also that the mean occupation num- 
ber does not considerably change over one oscillation pe- 
riod. The self-consistency of this approximation requires 
X 2 /{T ± cj r ) < 1. 

B. Results & discussion 

From Eq. (|2"0")) and Eq. (|2"Tj) we first of all see that up to 
second order in A our approach reproduces the LD cooling 
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rj | a | 17 1 a | 



FIG. 3: a) Cooling rate for a weakly driven TLS for different 
oscillation amplitudes a for Q = 0.05a»T- and T± = 0.15aj r . b) 
For the same parameters the cooling rate is plotted for fixed 
detuning as a function of a. The dashed lines indicate the 
analytic results given in Eq. (|24[) . c) The cooling rate r c (a) 
is plotted for different values of the Rabi frequency Q, S — 
— y/oJr ~ ^ 2 an d other parameters as above, d) Dependence 
of r c (a) on a for different values of T±. The values of Q and 
5 have been chosen to optimize r c (0). In all plots we have 
assumed Vu = A^tls = 0. 



rate and r c (a ->• 0) = T c (see Eq. (|A2]) in App. [A). 
The amplitude dependence of F c (a) is summarized in 
Fig. [3] for different parameters and assuming Fu, -/Vtls — > 
0. In Fig. [3] a) T c (a) is plotted for sideband resolved 
and weak driving conditions as a function of 5 and \a\. 
For resonant detuning S = — u r we find a degrading of the 
cooling rate with increasing \a\, which is expected from a 
saturation of the TLS. At larger amplitudes -q\a\ > 0.1 we 
also observe a strong modification of the cooling rates as 
a function of 5 and the appearance of additional peaks at 
5 = —2ui r , — 3uj r , . . . which correspond to multi-phonon 
transitions. For a better understanding of these features 
we change to the displaced oscillator basis introduced in 
Eq. ([7]) where the driving term in the Hamiltonian takes 
the form 

H n = - (a + e ? '( Qt - Q ) + (j_e- , '( Q+ - Q )) . (22) 

By expanding the exponentials we see that this Hamilto- 
nian couples different phonon number states \n) and for 
n ^> 1 the An-phonon transition matrix element is given 
by 

(l,n-An\Hn\0,n) =£ e -£ ^^W^AnOft 

ft _, 
— "2 JAn{2vVn)- 

(23) 



Here (x) denotes the generalized Laguerre polynomial 
and Jk{x) is the k-th order Bessel function. In the weak 
excitation limit £1 — > the cooling rate is proportional to 
the square of these matrix elements and for a detuning 
6 = — An x uj r we obtain a cooling rate 

r t ,a,^*A„x(£^)\ ,24, 

Fig. [3] b) shows that these analytic estimates correctly 
capture the dependence of T c (a) in the weak driving 
limit. 

In Sec. HI] we have seen that in the LD regime cool- 
ing rates are optimized under strong driving conditions 
Q, ~ uj r , and in Fig. [3]c) we study r c (a) for increasing 
Rabi frequencies. While the general dependence on the 
mean occupation number is similar to the weak driving 
limit, the cooling rates decrease much faster with increas- 
ing a and the benefit of strong driving is less significant 
at larger oscillation amplitudes. We attribute this ad- 
ditional suppression of the cooling rate to a deviation 
from the initial resonance condition, 5 = — y/cj^ — f2 2 , as 
the effective Rabi frequency changes with a. This effect 
is less severe for larger lincwidths which can be seen in 
Fig. El d) where we plot T c (a) for different TLS decay 
rates T±_ and with Q and 5 chosen to optimize r c (0). As 
we go from the sideband to the Doppler regime the cool- 
ing rates become more robust and the results from the 
LD regime are valid up to Xa < T±. Fig.[3]d) also shows 
that for Tj_ ~ u r and for certain values of a even small 
negative cooling rates can occur. 

In summary we find a significant reduction of the cool- 
ing rate already at moderately high occupation numbers 
T]yNth ~ 0.1 — 1. This degrading is most sever in the 
sideband resolved and strong driving regime, where in 
the LD limit cooling is optimized. However, for most pa- 
rameters and as long as S < we obtain r c (|a|) > 0,Va 
and an isolated resonator mode would still be gradually 
cooled towards the groundstate. Therefore, to under- 
stand the meaning of a temperature dependent cooling 
rate in a realistic setup we must study the full dynamics 
of the resonator mode and include environment induced 
re-thermalization processes. 



IV. FOKKER-PLANCK EQUATION 

In the previous section we have discussed the ampli- 
tude dependence of the cooling rate T c (a) under the as- 
sumption that the resonator mode is prepared in a coher- 
ent state | a). We now generalize this analysis to study 
the evolution of arbitrary resonator states p r (t). A natu- 
ral way to do so is in terms of a coherent state represen- 
tation for the density operator and in the present case 
a suitable choice is the Wigner function W(a,t) [39l ]. 
At high temperatures W(a = x + ip, t) corresponds to 
the classical probability distribution for position x and 
momentum p. In a frame rotating with u> r the Wigner 



7 



function evolves according to the Fokker-Planck equation 

(25) 

with a steady stately (a, oo) = e^ a]2 /{Nth+1/2) /{■n{N th + 
1/2)). To include the dynamics of the TLS we introduce 
three additional distributions 

W k (a,t) = ^J d 2 /3e- ia P*e- ia *PTr{e i(}at+i P* a a k p(t)}, 

(26) 

where k = —,+,z and {(Jk){t) = J d 2 aWk{a,t). The 
Wk(a,t) evolve on a timescale T\ and for A <C T±,uj r 
we can follow the arguments presented in Sec. IHII to de- 
rive an effective equation for the much slower dynamics 
of W(a,t). The details of this derivation is outlined in 
App. [B] and as a result we obtain the modified Fokker- 
Planck equation 

W(a, t) ~ i (^-a[iA(a) + T(a)} + H.c. ) W(a, t) 
2 \oa J 

( . . <9 2 M(a) a 2 M*(a) <9 2 \ T „. . 

(27) 

This equation describes the evolution of a damped res- 
onator where the damping rate T(a) = 7 + T c (a), the 
frequency shift A (a) and the diffusion terms D(a) and 
M(a) depend explicitly on the amplitude of the oscil- 
lation. The rate T c (a) is identical to the semiclassical 
cooling rate discussed above and is defined in Eq. (|2"l~|) . 
In the limit a — > we obtain for the diffusion terms 
D{a) = 7(7V tft + l/2)+r c (iVo + l/2) and M(a) = 0, such 
that Eq. (|2T|) correctly reproduces the cooling results in 
the LD limit. For general a the amplitude dependence 
of A (a), D(a) and M(a) can be calculated numerically 
as described in App. [Bj 

The parameters A(a), T(a) and D(a) depend on r = 
\a\ only while M(a) ~ a 2 M(\a\). Therefore, Eq. ([2"7]1 
preserves the radial symmetry of the initial thermal dis- 
tribution function and we can restrict the following dis- 
cussion to the radial equation, 

W(r,t) =1 (S~ r + 1 ) r (r)W(r,t) 

1 \ r 1 \ q (28) 

+** rr(p > (tf(r) + -) — W(r,*). 

Here we have made use of the fact that to lowest order in 
r\ the functions D(r) and M(r) commute with the radial 
derivatives and we have defined 

v ; r(r) 2 v y 

The steady state solution of Eq. (|28|) is 

"'(r.«-o.)-Af.-«W, BW;= y 

(30) 
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FIG. 4: Steady state occupation number n/ as a function of 
the thermal equilibrium value Nth- The parameters used for 
this plot are = 0.5cj r , 5 = — Vwp — fi 2 , r_i_ = 0.5cj r , 77 = 0.1 
and 7/wr = 0.25 x 10 _J . The value of N£ h indicates the 
crossover from the LD regime, rif = hld to a regime where 
cooling is strongly suppressed (shaded area) and n/ = JVt/,. 
The dotted line shows the approximate analytic expression 
for rif given in Eq. (|35[) . 



where M is a normalization constant and the final occu- 
pation number n/ is given by 

1 Z" 00 

nj = -- + 2nAT drr 3 e- R( - r) . (31) 

2 Jo 

A. Final occupation number 

In Fig. [4] we plot the characteristic dependence of the 
final resonator occupation number nj asa function of the 
initial value Nth- F° r very low Nth the final occupation 
number is limited by the intrinsic limit n/ = Nq and 
increases linearly n./ ps 7/r c (0) x iVt/, for moderate A^. 
This behavior is already captured by the LD result, nj ~ 
riic as discussed in Sec. |TTJ However, at very high N t h 
we see a crossover to a regime where cooling is no longer 
possible and tijkj Nth- 

To study this crossover in more detail we consider the 
relevant case where diffusive terms are dominated by 
thermal noise such that N(r)+l/2 ps yNth/T c (r). Indeed 
we find numerically that (D(r) + Re{M(r)}) < D(0), 
such that this assumptions is usually well justified, in 
particular for large r. Then 

R(r)^{^ + — f dr'r't c {r% (32) 
N t h n LD J 

where f c (r) = r c (r)/r c (0) and n LZ5 ~ 7 A/ t/l /r c (0) > 1 
is the final occupation number in the LD limit. From the 
discussion in Sec. lIIII we know that T c (r < 77) ~ r c (0) and 
in the limit Nth 3> n LD we obtain R{r < 77) ~ r^/ntn- 
In the opposite limit we can approximately write R(r ^> 
I/77) ps r 2 /N th +2i/ {r] 2 n LD ) where 

2 i = 2 ^ ctezl\Yr=^Y (33) 
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FIG. 5: Steady state Wigner-function showing the bimodal 
structure in the crossover regime Nth ~ Njr h . The parameters 
for this plot are Q/lj t = 0.1, S — — uv, T± — 0.15u; r , A = 0.1, 
■y/wr = 10 -5 and N t h = 2000. The dotted lines indicate the 
LD result and a thermal equilibrium distribution. 



is a numerical constant ~ Oil). Therefore, the steady 
state can be approximately written as the sum of two 
distributions, 



W{r) « N 



(34) 



as illustrated in Fig. [5j For ifriLD < T-\ the steady state 
distribution corresponds to the LD result while in the op- 
posite case r) tild > T-\ it approaches the thermal equi- 
librium. Therefore, while a degrading of cooling rates al- 
ready occurs at occupation numbers Nth ~ 1 /v 2 the final 
state of the resonator is still described by the LD theory 
as long as the much weaker condition n^u < Ii/rj 2 is 
satisfied. In particular, this implies that - within the 
validity of our model - ground state cooling is always 
correctly predicted by the LD approximation. 

From the approximate bimodal steady state distribu- 
tion given in Eq. (|3"4l we can calculate the final occupa- 
tion number 



information about the normalized cooling rate at a = 
as well as its temperature dependence. From Eq. (|37[) we 
see that N£ h has only a weak logarithmic dependence on 
the coupling strength A. 



B. Discussion 

For a weakly driven TLS we can use the analytic ex- 
pression given in Eq. (j2~4")l to show that X\ = 1 and 
I2 ~ £1 /oj%.- For strong driving SI ~ u> r and T± ~ uj r 
we obtain numerical values in the range of T 2 ~ 0.1 — 0.2 
which decrease again for T±_ -C ui r and T± 3> ui r . There- 
fore, assuming T± ~ uj r Eq. (|3"7) states that cooling of 
the mechanical mode with a TLS is only possible if the 
initial equilibrium occupation number is about 10-100 
times smaller than the mechanical quality factor Q. For 
cryogenic temperatures of T < 100 mK and mechanical 
frequencies of Lu r /2ir > 10 MHz we obtain Nth < 500 
and for reasonable Q-values of Q ~ 10 5 non-linear effects 
do not play a role. However, at T = 4 K or frequen- 
cies w r /27r < 1 MHz, the thermal occupation number 
Nth ~ 10 4 can be comparable with N^ h . Under these 
condition Eq. (|35l) must be used to determine the final 
occupation number for the specific set of experimental 
parameters. Finally, we see that at room temperature 
cooling will be difficult to achieve and requires the combi- 
nation of high mechanical frequencies with exceptionally 
high mechanical quality factors. 



V. MULTI-MODE EFFECTS 

Our analysis so far has been based on the model Hamil- 
tonian ^ where the mechanical resonator is approxi- 
mated by a single vibrational mode, e.g., the fundamen- 
tal bending mode. In reality the TSL is also coupled to 
higher order vibrational modes and more accurately 




1 - l/c 



LD 



X 1<LD 

1 + e N W IC LD , 



x Ni 



(35) 



where we have introduced the LD cooling factor Cld — 
r c (0)/7 3> 1. From this result we can deduce a critical 
thermal occupation number 



N th = S x 



c 



LD 



log(C 



LD 



(36) 



beyond which cooling is suppressed. In Eq. (|36[) N t 



th 



is written in terms of 77 and a given cooling factor Cld- 
However, both quantities scale like ~ A 2 and alternatively 



we can rewrite N^ h as 



H 



. (38) 

Here k = is the fundamental mode, a>o = a, which we 
have studied so far. The ak>o are operators for higher 
order mechanical modes of frequency uj^ which are cou- 
pled to the TLS with a strength A^. To provide a rough 
estimate of the influence of these higher order modes on 
the cooling rate of the fundamental mode let us switch to 
the displaced oscillator basis introduced in Eq. (J7J, but 
now with respect to all modes. We obtain 



n 
2 



(39) 



N, 



Ui 



r± log(Ciu) 



(37) 



where Q = tu r /^. Here we have introduced a single nu- 
merical constant I2 = (r c (0)rx/A 2 ) x X\ which contains 



where 77^ = Afc/w^ is the LD for the fc-th mode. For 
system parameters which are optimized to cool the fun- 
damental mode we can assume that non of the higher 
order modes is resonantly driven and therefore these 
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modes approximately remain in a thermal state. Then, 
by averaging Eq. Q39p over the thermal distribution of 
k > modes we obtain a reduced Rabi-frequency il —> 
f2(e-^ fc >° , ' fe ( a fc~ afc )) and a corresponding reduction of the 
zeroth mode cooling rate 

r c (a) -> r c (a)e-^ fc >o"^ 2Ar - +1 ) = r c (a) e -^*V 

Here we have expressed the reduction factor in terms 
of the LD parameter and equilibrium occupation num- 
ber of the fundamental mode and a numerical constant 
$ = 2 E fc >o(V A o) 2 x K°K') 3 - It depends on the ge- 
ometry of the mechanical beam [45[ and the exact cou- 
pling mechanism. For a TLS with a point-like coupling 
to the tip of a cantilever we obtain j3 f=s 0.0085, and 
(3 s» 0.011 for a TLS coupled to the center of a dou- 
bly clamped beam with low tensile stress. In both cases 
^>k ~ k 2 and only one or two high frequency modes con- 
tribute to the sum. For a high stress doubly clamped 
beam we obtain a slightly higher value of « 0.096, 
since the mode frequencies increase only linearly with k. 

The approximated cooling rate given in Eq. f|40[) shows 
that a single mode model for a mechanical beam coupled 
to a TLS is only valid for N t h < l/(f3r] 2 ). For higher tem- 
peratures the thermal motion of higher order vibrational 
modes has a significant influence on the system dynamics 
and in principle a full multi-mode model must be con- 
sidered. However, the influence of high frequency mode 
decreases quickly with increasing ui^ which suggest that, 
e.g., a two mode cooling scheme or an elaborate coupling 
design can be sufficient to eliminate these effects. 



VI. SUMMARY & CONCLUSIONS 

In this work we have analyzed the role of non-linear ef- 
fects on the cooling dynamics of a mechanical resonator 
coupled to a dissipative TLS. We have found that already 
at moderately high occupation numbers Nth ~ V a 
significant reduction of the cooling rate can occur. The 
reduction is most pronounced in the strong-driving and 
sideband resolved regime, where in the low temperature 
limit cooling is optimized. For the steady state we have 
found a crossover from the LD to the high temperature 
regime where cooling of a resonator mode with a single 
TLS is no longer possible. We have derived an approx- 
imate expression for the final occupation number which 
accurately describes this crossover and discussed the de- 
pendence of the critical occupation number N£ h on the 
relevant system parameters. Our results are relevant for 
experimental implementations of cooling schemes oper- 
ating at temperatures above ~ IK and in general for 
cooling of resonator modes with frequencies of ~ 1 MHz 
and below. 

The theoretical approach presented in this paper pro- 
vides a unified description of the dynamics of a weakly 
coupled resonator-TLS system in the quantum as well as 
the high temperature regime. It is valid for A <C Tj_,u r 



which is fulfilled for many solid state TLS coupled to 
nano-mechanical resonators. Therefore, apart from cool- 
ing our analysis can be adopted to study related prob- 
lems like non-linear dissipation mechanisms caused by 
single or few two-level fiuctuators [46[ or phonon lasing 
physics Ha]. 
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Appendix A: The fluctuation spectrum S(u) 



To evaluate the fluctuation spectrum S(u) defined in 
Eq. dTTJ) we write S(u) = A 2 /2 x Re{C 3 (s = -iuj)} where 
the vector C(s) is the Laplace transform of the set of 
correlation functions C(t) = (S(t)(Tz)q — (5)o(f a )o and 
S = (<7_, <t + , a z ) T . Using the quantum regression theo- 
rem 39] we obtain 



C(s) 



1 



si 



(Oo 

-(or+)o | „ <t, „ 



(Al) 



where the matrix A is defined in Eq. (fT8|) and the steady 
state expectation values (S)q follow from the Bloch equa- 
tions of the TLS, (S) = A (S) - r ± U z , where V z = 
(0, 0, 1) T . The cooling rate can be further simplified and 
written as 



r c = -a 1 



which agrees with the result derived in Sec. IIIII for 
T c (a -► 0). 



Appendix B: Fokker-Planck equation 

In this appendix we derive the effective Fokker-Planck 
equation (1271) . We change into a frame rotating with 
the mechanical frequency uj r such that the free evolution 
of the resonator Wigner-function W{a,t) is W(a,t) — 
T> y W(a,t). The differential operator for mechanical 
damping 2? 7 is defined by Eq. (1251) . The distributions 
Wk(a,t) introduced in Eq. (|2"6"|) are grouped into a vec- 
tor W = (W-, W+, W Z ) T which for A — > evolves as 

VV(a, t) = AW(a, t) - T ± V z W(a, t) + V 7 W(a, t). (Bl) 
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For finite A the evolution of W(a) and W(a) is coupled 
by the TLS-resonator interaction H\(t) = ■|(ae _i " r * + 
a) e lUrt )a Z: and we obtain the additional terms 

W ah (a) = - i± {ae^ + a*e^) W Wk ,^(a) 

_ ( e -i"rtJL _ j"rtji\ Ws . ( a ) 

■ 1 Q a * e ) ^{<r*><M+W' 



<9a 



(B2) 

where {., .}+ is the anti-commutator. To remove the ex- 
plicit time dependence we change to Floquet representa- 
tion, 

oo 

W k (a,t)= £ Wj!{a,t)e"^ nt , (B3) 

n— — oo 

and obtain the set of coupled equations 

W n (a, t) = iuj r nW n {a, t) + V^W n (a, t) 

+4(^" + v<>-^»'rW)), (B4) 

and 

W>, t) = A n W n (a, t) - T ± V Z W n (a, t) + V 7 W n (a, t) 
- i\aA z W n - x (a,t) - iXa* A z W n+1 (a, t) 



d 



d 



—W n -\a,t) - —W n+1 (a,t) 



da 



da 



(B5) 

where we have defined A n := A + iuj r n\. Our goal is to 
eliminate the fast dynamics of the TLS and to derive an 
effective Fokker-Planck equation for the slowly varying 
Wigner- function W°(a,t). To do so we point out that 



X^W k (a,t) 



Ac 



n r (t) 



W k (a,t) ~ 0(X)xW k (a,t), (B6) 



where n r (t) is the mean resonator occupation number. 
Therefore, in the regime A <C r_i_,u; r we can formally 
expand the coupled differential equations in powers of 
A x d/da. 

Cooling rates. To zeroth order in A x d/da we ob- 
tain for the resonator mode W°(a,t) = T> 7 W°(a,t) and 
W n (a,t) — otherwise. In Eq. (IB5|) we must include 
terms ~ a A and to zeroth order in our expansion pa- 
rameter we obtain W£(a, t) = S£ (a) x W°(a, t). Here 
the SJ^ (a) have been introduced in Sec. IIIII and they 
are defined by the matrix continued fraction in Eq. (1201) . 
Note that in deriving this expression we have neglected 
the influence of mechanical damping. If the resonator 
state is close to thermal equilibrium T> 7 Wk ~ 0(7), 
while far away from equilibrium T> 7 Wk ~ 0(N t h^/)- 
However, to obtain a non-equilibirum state we require 
j Nth <S r c < T± and therefore as long as T± ^> 7 this 
approximation is automatically justified. 

After inserting the zeroth order results for Wu(a) back 
into Eq. (|B4[) we obtain to first order in A x d/da, 



W°(a,t) = V 7 W°(a,t) 



l 2 



W°(a,t). 



(B7) 



From this expression we can deduce a cooling rate F c (a) 
and a frequency shift A (a) by setting 



i\S+l(a) =: a (T c (a) + iA{a)) , 



(B8) 



which reproduces the semiclassical results derived in 

sec. inn 

Diffusion. We now include corrections to W£ (a) which 
are first order in A x d/da and lead to diffusion terms in 
the equation for W°. This requires some care and we 
first point out that apart from damping the zeroth order 
results also leads to a static force ~ S z0 (a). In the LD 
regime S® (a) = (o~ z )o and in Sec.|H]we have eliminated 
this force by redefining the resonator equilibrium posi- 
tion. To obtain physically meaningful result the same 
must be done here and therefore we perform a variable 
transformation a — > a + X/(2u> r )S z Q (a)e lUrt . To lowest 
order in A this transformation does not change the result 
given in Eq. (|B7|) . However, in Eq. ()B5|) this transforma- 
tion introduces two additional terms 



ri(a,t) = ... + i±S%(a)-^W n - 1 (a,t) 
-i^S° Zt0 (a)-^W n+1 (a,t), 



(B9) 



which must be taken into account to correctly reproduce 
the diffusion terms found in the LD regime. 

We can now insert the zeroth order result Wu(a,t) = 
S£ (a) x W°(a,t) back into the modified Eq. (|B"5j) 
and iteratively solve for corrections which are linear in 
A x d/da. Note that ^S n (a) = S n {a)-J^+(D(X) and to 
lowest order in A we can exchange S^(a) and the deriva- 
tives. Then, the final result can be written in the form 

iAW+Vi*) = a(r c (a) +iA(a))W°(a,t) 

d d (BIO) 

+ f2(a) — W (a,t) + f 3 (a)-^W (a,t). 



After inserting this expression back into Eq. (|B4|) we end 
up with the Fokker-Planck equation (|27p. where D(a) := 
Re{f 2 (a)} + j{N th + l/2)&ndM(a) := f 3 (a). In the LD 
regime, i.e. up to second order in A, we obtain 



i\W+ l {a,t) = -X 2 a(0,0,l)(A +1 )- 1 A z S W°(a,t) 

(Bll) 



y(0,0,l)(A +1 )- 1 (( ( 7 z ) S ? o-K) JLw°(a,t). 



By comparing this expression with the result derived in 
App. \M we see that r c (a) = T c and D(a) = T C (N + 
1/2) + j(N t h + 1/2) agree with the LD results. For ar- 
bitrary a the coupled set of Eqs. (|B5[) can be evaluated 
numerically by truncating the coupled set of equations at 
a large value of \n\. 
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